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^ Abstract 

^ The recently discovered iron-based superconductors AyFe2-xSe2 {A=K, Rb, Cs, Tl) show a long- 

range antiferromagnetic order with an unexpected high transition temperature Tjv ~ 550 K and a 
large effective magnetic moment ~ 3.3/2^ at iron sites. This unusual behavior is arguably related 
J, ■ to the existence of a unique vacancy order. Taking the extended J1-J2 Ising spin lattice 

^ . as a minimal model, we investigate the finite-temperature magnetic phase transitions in a square 

lattice with a \/5 x -v/5 vacancy superstructure by using large-scale Monte Carlo simulations. 
By the parallel tempering technique, the block spin checkerboard and stripe antiferromagnetic 
states are detected to be the groundstates for two simplified sets of model parameters. The short- 
time dynamic approach is utilized to accurately determine the critical temperature as well as the 
static and dynamic exponents. Our results demonstrate a dramatic enhancement of the critical 
temperature in the non-frustration case due to the vacancy order. 
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I. INTRODUCTION 



The discovery of superconductivity in iron pnictides [l-5| has renewed an intensive study 
on the interplay between superconductivity and antiferromagnetism {g]. A broad family of 
the iron-based superconductors has been synthesized, by which the parent coinpounds are 
typically represented by the 1111-type LaFeAsO l|, the 122-type BaFe2As2 j?], the Ill- 
type LiFeAs [8], and the 11-type FeSe [9]. The parent compounds show antiferromagnetic 
transitions around the Neel temperature T/v ~ 100 — 200 K. The long-range magnetic order 
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is a stripe-like (or collinear) antiferromagnetic state 
magnetic order is a bi-collinear antiferromagnetic state 
as superconductivity are closely related to a common two-dimensional (2D) Fe-atom square 



, ll| except for the 11-type where the 
[l2[. The magnetic properties as well 



lattice 



13Nl5|. The effective magnetic moments of each irons determined by experiments 
e— l.O/is/Fe, while the iron moments estimated from the first- 



are usually within 0.3/xr/ 



principles calculations 



16 



- |l9j or model analysis 20| could be around 2.0/iB/Fe or larger. 



Recently, a new family of iron-based superconductors, i.e., the intercalated iron chalco- 



genides AyFe2-xSe2 {A=K, Rb, Cs, 
transition temperature T^c ~ 30 K 



Vl has been found with moderate high superconducting 
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23| . These materials are structurally similar to the 



122-type iron pnictides, thus leaving certain amounts of Fe- vacancies in the iron square sub- 
lattice. The iron vacancies are expected to order in some periodic superstructures, rather 
than to distribute randomly within the FeSe layer, resulting in the normal state insulating 
behavior 22|. Indeed, the a/5 x \/5 vacancy ordering pattern as shown in Fig. [1], correspond- 
ing to X = 0.4, seems to be most stable as confirmed by the neutron diffraction 2J-|26| and 
transmission electron microscopy 27| experiments. In addition to the a/5 x a/5 Fe-vacancy 
superstructure, a novel magnetic ordering pattern, i^ block spin checkerboard (BSC) state 
is observed by the neutron diffraction experiment 2J]. Moreover, the magnetic transition 
temperature is unprecedentedl y h igh, reaching T^r ~ 5597^^, and the effective magnetic mo- 
ment is as large as 3.31/iB/Fe \24\. 

The early theoretical explorations to the magnetic and electronic structures of AyFei,QSe2 
based on the first-principles calculations have revealed the BSC state as the groundstate with 
effective iron moment ~ 2.8/iB/Fe— 3.4/i5/Fe, and a band gap ~ 500 meV a.ty = 0.8 



28 



Another silent feature revealed by these calculations is a significant lattice contraction of 
the fundamental iron blocks but without breaking the symmetry of lattice structure. The 
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observation leads to a microscopic consideration for the magnetic structure based on the 
extended J1-J2 spin model |28l | , -w here the BSC state could be the groundstate for certain 
range of model parameters 28|, l30| • This model involves the nearest-neighbor (NN) and the 



next-nearest-neighbor (NNN) exchange interactions between iron moments, and captures the 
vacancy superstructure and the lattice distortion in a minimal manner. Other models like the 
J1-J2-J3 spin model, which emphases the relaxation of magnetic frustration by the third- 



31| or with 32] a biquadratic interaction 



nearest-neighbor exchange interaction without 
term, can also account for the BSC state. 

Experimentally, the vacancy order and the BSC state coexist with the superconductivity 
in the AyFe2-x'^^2 compounds for y > 0.8, x ^ 0.4. This raises heated debate on whether 



the co-existence is _an intrinsic property of a single phase electronic structure 
phase separation 



33| or due to a 



34| . Theoretically, the importance of the vacancy order and the BSC state 



in the 
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brmation of the insulating phase as well as superconductivity has been investigated 



42| . However, the question that why the magnetic transition temperature T^r is 



much higher than the ones in other iron pnictides or chalcogenides, has not been addressed 
so far by detailed calculations. 

This fundamental question, though seems obvious at first glance, is actually non-trivial. 
Firstly, it is believed that unconventional superconductivity, which frequently emerges near 
the border of antiferromagnetic phases such as in the cuprates, heavy fermions, and iron 
pnictides/chalcogenides, is usually related to the magnetic fluctuations with a characteristic 
energy scale roughly proportional to the magnetic transition temperature Ttv 43|. But the 
present class of intercalated iron chalcogenides is an obvious exception. Secondly, while 
the intralayer magnetic interactions Ji and J2 in AyFe2-3^2 is comparable to those in 



iron pnictides 



44| as well as in 74Fe2Se2 without vacancies 45|, |46| , the interlayer magnetic 



coupling is very small in AyFe2-x^^2 |28[ and is not sufficient to account for the enhancement 
of T/v . Thirdly, it is well-known that the vacancies in spin systems always lead to a reduction 
of transition temperatures such as in the randomly site-diluted Ising model j47]. Therefore, 
detailed calculations are required to understand why and how the opposite tendency could 
appear in the AyFe2-x'^Q2 compounds. 

In this paper, we consider this question by performing systematic simulations on the 
finite-temperature magnetic phase transitions in the extended J1-J2 spin model defined in 
a square lattice with a -\/5 x y/h vacancy superstructure. This model has a rich phase 
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diagram consisting of numerous groundstates when tuning the model parameters 28|, |30 |. 
We then focus on three simphfied yet representative sets of model parameters. Our goals are 
of two-folds: to clarify the nature of the finite-temperature magnetic order- disorder phase 
transitions and to understand the mechanism driving the transition temperatures. 

There are several technique difficulties to our numerical investigations. For the systems 
with ( Ji, J2) couplings and vacancies, sufficiently large system sizes are required in order 
to find the true groundstate among numerous possible magnetic configurations. For this 
purpose, we find that the parallel temperin g te chnique based on Monte Carlo method, which 



was applied in studying spin glass systems! 



is an appropriate approach to searching for 



the groundstate of the present model. Once the groundstate is determined, the magnetic 
order- disorder phase transition at finite temperatures can be investigated by large-scale 
Monte Carlo simulations developed for the systems with vacancy superstructures. Due to 
severe critical slowing down, however, accurate determination of the transition temperature 



and critical exponents is very hard. In this respect, the short-time dynamic approach 49 



can be utilized. It is particular efficient if we start from the lower temperature ordered phase. 
Recent activities include various applications and developments 5lL l52| such as theoretical 



and numerical studies of the Josephson-j unction arrays 



53| and ageing phenomena 
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55|. 



Very recently, the depining transition and the relaxation-to-creep transition in the domain 
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-|59|. 



wall motion have been investigated within this approach 

With the parallel tempering technique and the short-time dynamic approach, we are able 
to calculate the groundstate as well as to accurately determine the transition temperature 
and critical exponents. In Sec. II, the model and scaling analysis are described, and in Sec. 
Ill, the numerical results are presented. Finally, Sec. IV is devoted to the conclusions. 



II. MODEL AND SCALING ANALYSIS 



A. Model and methods 



The extended J1-J2 model is the extension of standard J1-J2 model defined on the Fe- 
square lattice of iron pnictides {20! to the situation with a Fe-vacancy superstructure. Here, 
the distance between two NN vacancies is in unit of the distance between two NN 
iron sites. Owing to the -\/5 x a/5 vacancy superstructure, the whole lattice consists of 
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the fundamental minimal square blocks containing four iron sites. Due to the symm etry 
invariant lattice distortion, the intrablock and interblock interactions could be different 28 |. 
The model is then defined by the following Hamiltonian, 

H = '^^{JlSn,aSn,a+l + J2Sn,aSn,a+2) 
n.a 

•n,a n,a 
+ Sn,asSn-\-S~l,as-l)i (l) 

where n denotes the block index, n + S is short for the nearest-neighboring block, a is the site 
index which goes from 1 to 4, and as selects the site connecting to the nearest-neighboring 
block. Jl and J[ ( J2 and Jj) are the NN (NNN) couplings of intra- and interblock, as shown 



in Fig. [H Since an almost saturated magnetic moment at the iron site is reported |24j . 
the quantum effect is suppressed at finite temperatures due to the large local spin. So the 
corrections due to quantum fluctuations to the magnetic properties at finite temperatures 
can be safely neglected. Therefore, the spin is treated as a classical Ising spin, i.e., Sn,a — il. 
The specific value of IS"! is not crucial in the following discussions. 

The general phase diagram of this classic J1-J2 spin model is very complicated. A previous 
Monte Carlo study suggests a phase diagram with several specified groundstates for the 
vector version of this model with certain model parameters [30]. For our purpose, we consider 
more simplified cases by fixing the non-zero coupling strengths \J\ = 1. Explicitly, three 
sets of model parameters are considered, 

Case-I: J^ = J^ = -1, J[ = J'^ = l 

Case-II: = J2 = 1, J[ = J'^ = l . (2) 

Case-Ill: = -1, J2 = 0, J( = 1, = 

Here, Case-I involves the ferromagnetic intrablock NN and NNN couplings and antiferro- 
magnetic interblock NN and NNN couplings, while all NN and NNN couplings in Case-II are 
antiferromagnetic. Case-Ill is a simplification of Case-I, i.e., without the NNN interactions. 
First-principles calculations suggest that in realistic AyFe2-x'^Q2 compounds, Ji, J2, and J[ 
are ferromagnetic, while is antiferromagnetic. However, J!^ dominates over J[ not only 
because Jn is larger than I J,' I , but also the number of interblock NNN sites are twice of the 

n 

interblock NN sites |2g]. As we will see later, the sign of J[ is not crucial in this case. The 
origin of the ferromagnetic couplings is due to a combined effect of Hund's rule coupling and 



short-ranged hopping integrals (of Fe 3(i-orbitals and Se 4p-orbitals) which are enhanced by 
the vacancy-induced lattice contraction. Notice that when no iron vacancies appear as such 
in AFe2Se2 
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46| . both NN and NNN interactions are antiferromagnetic as in other iron 



pnictides. 

In the following, we shall find the magnetic configurations of the groundstate and inves- 
tigate the finite-temperature magnetic order-disorder transitions in each of the three cases 
with Monte Carlo simulations. To overcome the critical slowing down around the phase 
transition, we adopt the short-time dynamic approach. Usually, two relaxation processes, 
i.e., those starting from the groundstate (ordered state) and high temperature state (disor- 
dered state) are considered. To extract the transition temperature and critical exponents, it 
is more efficient to study the dynamic relaxation starting from the groundstate. Fig. [T] shows 
two typical magnetic configurations, corresponding to the magnetic ordering patterns of the 
BSC state and the stripe (or coUinear) antiferromagnetic state (denoted by SAFM). The 
two magnetic ordering patterns are respectively observed in realistic AyFe2-xSe2 compounds 



and iron pnictides by experiments 10 



24| . We will show that the groundstates of Case-I 



and Case-II are the BSC and SAFM states, respectively. 

For a simple model, such as Ising model with only NN couplings on a square lattice, it is 
straightforward to obtain the groundstate magnetic configuration according to the symme- 
try. However, the vacancy order and frustrated antiferromagnetic NNN interactions in the 
extended J1-J2 model make this task difficult, because the number of competing magnetic 
configurations may increase rapidly with the system size. Large-scale simulations imple- 
mented by the parallel tempering algorithm are then performed to find the true groundstate. 
The details of the parallel tempering algorithm can be found in Ref. 48|, and the main idea 
is briefiy illustrated below. In this algorithm, m parallel replicas are analyzed, each of which 
is performed independently at a fixed temperature Tj (Ti < Tj < T^). Following the refer- 
ence, we fix Ti = 0.1, Tm = 1.6 and set T^+i — Tj = {T„i — Ti)/ (m — 1). In order to avoid the 
situation where replicas at low temperatures get stuck in local minima, one can swap the 
configurations of two randomly selected temperatures Tj and Tji. Starting from a random 
initial condition, a standard Monte Carlo dynamics is performed in each replica, and a trial 
exchange of two configurations Xj and Xji (corresponding to the jth and j'th replicas) is 
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attempted periodically, and accepted with the probability 

I exp(-A), for A > 
W{X„K,\X,,,K,,) = \ . (3) 

[ 1, for A < 

Where A = —{Kj — K'-){Hj — H'-) is defined with the inverse temperature Kj = l/Tj 
and Hamiltonian energy Hj. For convention, we restrict the replica exchange to the case 
j' = j + As time evolves, the magnetic configuration at the lowest temperature approaches 
to the groundstate. 

After preparing the groundstate as the initial state, we update the spins with the heat- 
bath algorithm. Our simulations are performed with lattice sizes L = 250, 500, 1000, up to 
tmax = 25,600 Monte Carlo step (MCS). Here MCS is defined by L x L single-spin flips. 
Periodic boundary conditions are used along the x and y directions, respectively. For each 
case, more than 16, 000 samples are performed for average. Errors are estimated by dividing 
the samples into three or four subgroups. If the fluctuation of the curve in the time direction 
is comparable with or larger than the statistical error, it will be taken into account. 

To investigate the dynamic relaxation, the pseudo-magnetization M{t) = M^^^t) and its 
second moment M*^^)(t) are introduced by the projection to the groundstate. 



A; = 1,2, (4) 



where Si{t) is the spin at the time t on the lattice site i, Xi is the one from the groundstate, 
L is the lattice size, and (■ • •) represents the statistical average. The pseudo-magnetization 
M{t) plays a role as the order parameter of the magnetic transition. When the groundstate 
is degenerate, a computationally convenient root-mean-square order parameter is introduced 
6o| . Other important observables are the susceptibility and Binder cumulant U{t), 

x{t) ~ M(2)(^)-M(^)^ 

U{t) ~ x{t)/M{tf. (5) 

For the dynamic relaxation starting from the disordered state, the spatial correlation func- 
tion C{r,t) and two-time correlation function A{t,t') are measured, 

AM = E (5.(05. W), (6) 

i 

where r is the spatial distance, t' is the waiting time, and ci = 2 is the spatial dimension. 



B. Scaling analysis 



The magnetic order- disorder transition at finite temperatures in the present model is of 
;he second order, compatible with the magnetic transitions in the 74j^Fe2_a;Se2 compounds 



24j |. where no lattice structural transition is accompanied with the magnetic ordering except 
for the iron vacancy ordering stability taking place at an elevated temperature Ty ~ 580 K. 
Hence, one expects that the order parameter M{t) should obey the dynamic scaling form, 
after a microscopic time scale tmic l^sj], 

M^''\t,T,L) = t-''^/'''M{t^/^'r,t^/'/L), (7) 

here (3 and u are the static exponents, z is the dynamic exponent, and r = (T — Tc)/Tc is 
the reduced temperature. denotes the transition temperature which can be either Curie 
temperature Tc in the ferromagnetic transition or Neel temperature T/v in the antiferromag- 
netic transition. On the right side of the equation, the overall factors t~^^/'^^ indicates the 
scaling dimension of M{t), and the scaling function M{t^^'^^T, t^^^ /L) represents the scale in- 
variance of the dynamic system. For a sufficiently large lattice and in the short-time regime, 
the nonequilibrium spatial correlation length ^{t) ~ t^/^ is much smaller than the lattice 
size L. Therefore, the finite-size effect is negligible, and a power law behavior is expected 
at r = 0, 

M(t) ~ r^'^^ (8) 

With Eq. ([7]), the precise location of the transition temperature Tc is determined by searching 
for the best power-law behavior of M{t, r), and the critical exponent 1/vz is measured from 
the time derivative of lnM(t, r). 

For the susceptibility, it is different. Since ^{t) is small, the spatially correlating terms 
(5*1 Xi 5*2X2) with |ri— r2| > ^{t) can be neglected. In other words, one of the two summations 
over ri and r2 in M^'^\t) — M{tY is suppressed. It then leads to the finite-size behaviors 
xit) ~ L^'^. Together with Eqs. ([7]) and ([8]), one may derive the scaling forms, 

Xit) ~ t^'^^/L\ 

U(t) ~ t'^'^/L'^, (9) 
with the scaling law --y/u = d — 2j3/u. 
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For the dynamic relaxation starting from the disordered state, the correlation functions 
C{r,t) and A{t,t') should obey 

A{t,f) ~ t'-'^/'^'Am/m), (10) 

where C{s) and A{q) are the scaling functions with s = r/^{t) and q = C{t)/C{t')- Together 
with Eqs. and (fTUj) . one may derive the scaling form of the integral S{t) = J C{r,t)dr, 

S{t) ~ t('io-2/3/.)/.^ 

where do denotes the dimension of the integration. For a sufficiently large lattice and at 
the critical point Tc, a surprising increasing behavior of the pseudo-magnetization M{t) is 
observed, 

M{t) ~ niof. (12) 

Here mo is the initial rnagnetization, and ^ is a local critical exponent, reflecting the effect 
of the initial condition [61 1. 



III. MONTE CARLO SIMULATIONS 



As shown in Fig. [2]^a), the second moment of pseudo- magnetization is displayed in the 
parallel tempering process at the lowest-temperature replica for Case-I and Case-II. The 
lattice size L = 200 and replica number m = 10 are used. As time grows, the curves 
approach to the unit. It indicates that the BSC and SAFM states, displayed in Fig. [H are 
the true groundstates for Case-I and Case-II, respectively. As a test, we also consider a set of 
model parameters by re-scaling the couplings obtained from the flrst-principles calculations 
In this case, more than 1000 samples are performed, up to t^ax = 1? 000, 000 MCS. All 



of them evolve to the BSC state. This reflects the fact that there is an extended region of 



the BSC state in the groundstate phase diagram of the model 



30|. 



A. Dynamic relaxation from groundstate 

With Monte Carlo simulations, the dynamic relaxation starting from the groundstate is 
investigated. In Fig. Mjo), the time evolution of the pseudo-magnetization M(t) in Case-I 

9 



is displayed for different inverse temperatures K = 1/T with the lattice size L = 500. The 
curve drops rapidly down for smaller K, while approaches a constant for larger K. Searching 
for the best power-law behavior, the critical point Kc = 0.27595(3) is determined accurately. 
According to Eq. ([H]), one measures the exponent P/uz = 0.0585(6) from the slope of the 
curve at Kc- Additional simulations with L = 250 and L = 1000 confirm that the finite-size 
effect is already negligibly small. For comparison, the dynamic behavior in Case-II is also 
studied, and the critical point Kc = 0.8148(1), the exponent (3/h'z = 0.0556(3) are derived. 

In order to approximate the differentiation of In M{t, r), the simulations at temperatures 
in the vicinity of the critical point are performed. In Fig. ^a), a power- law behavior of 
the curves is observed but with certain corrections to scaling at the early times. A direct 
measurement from the slope gives the exponents 0.471(5) and 0.512(3) for Case-I and II, 
resp ectively. After introducing a power-law correction to scaling, dr InM (t) ~ t^^'^^il + c/t) 
|55| . one can fit the numerical data extending to rather early times. It yields l/uz = 0.468 
in Case-I, and 0.510 in Case-II. 

In Fig. [3t^b), the time evolution of the Binder cumulant U{t) is plotted at for the 
two cases. The possible finite-size behavior is also investigated with different lattice sizes 
L = 250, 500 and 1000, and data collapse is observed according to Eq. (|9]). From the slope, 
one measures the exponent d/z = 0.928(5) in Case-I, and 0.921(5) in Case-II. 

Finally, according to the measurements of P/i^z, l/i^z, and d/z, we calculate the indi- 
vidual exponents /3 = 0.125(2), = 1.00(1), z = 2.16(1) in Case-I, and /3 = 0.109(1), z/ = 
0.90(1), z = 2.17(1) in Case-II. 

B. Dynamic relaxation from disordered state 

Now we turn to the dynamic relaxation starting from the disordered state at the critical 
temperature Tc. In Fig. Hl^a), the spatial correlation function C{r,t) is displayed for Case-I 
as a function of distance r at different time t. To confirm the scaling behavior of C(r, t), for 
example, we fix t' = 20480 MCS, and rescale r to {t'/ty/'r and C{r,t) to {f /t)-^l^^^'C{r,t). 
Data of different t nicely collapse to the curve of f with the exponents P/uz = 0.0585 
and z = 2.16 as input. A power-law decay is then observed at small s = r/^{t) with the 
slope 2/3/z/ = 0.25(1). In order to extract the characteristic of the scaling function, C(s)s°'^^ 
against s is plotted in the inset. For large s (e.g. s > 2), an exponential behavior is detected. 
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indicating the scaling form, 

C{s) ~ 5-2/5/- exp(-as). (13) 

Together with Eqs. ( ITOj) and (fT3|) . one may derive the critical behavior of the spatial corre- 
lation function C{r,t) in the limit r/C,{t) — ?■ oo, 

Cir,t)r^-^expi-ar/at)), (14) 

here ^(t) ~ t^/^ is the spatial correlation length. 

In Fig. m^b), the integrated correlation function S(t) is displayed for Case-I, and the 
exponent {do — 2f3/h')/z = 0.331(8) is estimated from the slope, according to Eq. (ITT]) . The 
dimension do = 0.97(2) is calculated, very close to 1. Similarly, a power law behavior is also 
observed for the susceptibility x(t) with the slope -y/uz = 0.806(4). It yields the exponent 
7 = 1.72(2). While similar measurement for Case-II yields a different value of 7 = 1.58(2). 

The scaling behavior of the two-time correlation function represents a kind of ageing 
phenomena MjIsSj. According to Eq. (ITU]) , the scaling function A{t/t') is plotted in Fig.[n](a), 
as a function of g = ^{t)/^{t'). Obviously, data for different waiting time t' collapse onto a 
master curve, and exhibits a power law decay in the large q regime (e.g. q >2). It indicates 
that the scaling function A{q) takes the form 

A{q) ~ q-\ (15) 

with the scaling law A = d — 9z. According to the formula, the critical exponents A = 1.59(1) 
and 1.65(1) are estimated for Case-I and II, respectively. 

Finally, a surprising increase of the pseudo-magnetization M{t) is displayed in Fig. El^b) 
with the lattice size L = 1000. From the slope of the curve, one measures the critical expo- 
nent 6. Strictly speaking, 6 is defined at the limit mo — ?■ 0. However, practical measurements 
at this limit is not possible. In this work, the initial magnetization mo = 0.01 is prepared, 
which is believed to be small enough. It yields the exponent 6 = 0.186(2) in Case-I, larger 
than the corresponding value 0.167(1) in Case-II. 



C. Discussion 



All the measurements of the transition temperature and critical exponents are summa- 
rized in Table [H in comparison with those of the 2D square Ising model without vacancies. 
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In general, the vacancies would lead to a continuous decrease in the transition tempera- 
ture |47|. Magnetic frustration induced by the antiferromagnetic NNN interaction will also 
decrease the transition temperature 60|]. For example, the calculated critical temperature 
T(. = 1/ Kc. = 1.2273(1) in Case-II is much lower than the one 2.2692 in the 2D Ising model. 
The reduction of should be due to both vacancies and magnetic frustration. However, 
a dramatic enhancement in the transition temperature, Tc = 3.6238(4), is found in Case-I. 
This value is much larger than that of the 2D Ising model, and almost three times as large 
as that in Case-II. Differences between Case-I and Case-II are also observed in individual 
critical exponents /3, 7, 6, and A which differ by about 10%. It suggests that the magnetic 
transitions in Case-I and Case-II are not in the same universality class. Further compari- 
son shows that the former belongs to the Ising universality class, while the latter does not. 
Interestingly, the ratio /3/z/, 7/z/ and the dynamic exponent z in Case-II agree well with the 
corresponding values of the 2D Ising model. It supports the dynamic generalization of the 
"weak universality" hypothesis proposed by Suzuki [60,162], where only the reduced critical 
exponents 7/z/ and z = /\/v are universal, irrelevant to the details of the interactions. 

In order to understand above results, a simpler example, Case-Ill defined in Eq. ([2]), is 
investigated. Using the parallel tempering algorithm, the BSC state is confirmed as the 
groundstate, too. In Fig. [6]^a), the inverse transition temperature = 0.6952(1) and the 
exponent P/i^z = 0.0564(4) are measured from the dynamic relaxation from the groundstate. 
Other critical exponents are also calculated, as shown in Table. [B Except for the transition 
temperature, the critical exponents in both Case-I and Case-Ill are very close to the ones 
of the 2D Ising model, showing that they are both in the Ising universality. 

Now we argue that the agreement of critical exponents for Case-I, Case-Ill and the 2D 
Ising model is not an accident. The present model, though with vacancies, has a perfect 
symmetry that each site has three equivalent neighbors preserving both -\/5 x trans- 
lational and four- fold rotational invariances of the lattice structure 4l|. In particular, 
Case-Ill is invariant under a block-spin rotation: Si — SiC^^'", associated with a mapping 
Ji — )■ Ji, J[ — )■ — J(. Here, the integer n denotes the block index, and Ji (J() indicates 
the coupling of intrablock (interblock). Therefore, the model is equivalent to the ferromag- 
netic Ising model defined on the square lattice with the vacancy superstructure. Then, a 
topological deformation can be performed from the square lattice with the a/5 x ^/S va- 
cancy order to the bathroom-tile lattice, as illustrated in Fig. El^b). Hence Case-Ill and the 
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2D bathroom-tile ferromagnetic Ising model are equivalent. Remarkably, the latter model 
(with the NN coupling only) is exactly solvable, with the exact inverse Curie temperature 



= tanh-^{^(b + 4v/2)/2 - 1 - 1/^2) ^ 0.6951 [63|. This value is in perfect agreement 
with our numerical value 0.6952(1). In addition, it is known that the Ising models defined 
on the square, triangular, Kagome, and bathroom-tile lattices belong to the same universal- 



ity class j49|, l64j-l67j. As a consequence, identical critical exponents are predicted between 



Case-Ill and the 2D Ising model, as revealed in our numerical results. We note that the 
critical temperature T^. is lower in Case-Ill than in the 2D Ising mode. This is clearly due 
to the existence of 20% vacancies (x = 0.4 in our model) which in turn leads to three NN 
bonds for each iron site. 

Similar analysis can also be carried out for Case-I and Case-II. Under the block-spin 
rotation and the topological transformation, the bathroom-tile Ising model with both the 
ferromagnetic NN and NNN couplings is derived for Case-I. Since the ferromagnetic NNN 
coupling is irrelevant to the Ising universality, the critical exponents agree well with those of 
the 2D Ising model. However, it significantly increases the transition temperature due to 
the increase of the ferromagnetic coupled bonds. By contrast, the situation is quite different 
in Case-II, because the magnetic frustration between the antiferromagnetic NN and NNN 
couplings exists still even after the mapping. It explains why the individual exponents are 
non-universal (may vary with model parameters in the same SAFM phase) and different to 
those of the 2D Ising model. Meanwhile, the transition temperature is also suppressed. 

We note that for the parameters obtained from the first-principles calculations, i.e., Ji, J2, 
and J[ are ferromagnetic, is antiferromagnetic, above mapping leads to dominating ferro- 
magnetic J 2 and small antiferromagnetic J[. Therefore, the magnetic frustration is actually 
suppressed by dominating J'^, leading to the same BSC state as in Case-I. 



IV. CONCLUSION 



Using the parallel tempering technique and short-time dynamic approach, we have numer- 
ically investigated the finite-temperature magnetic order- disorder transitions in the extended 
J1-J2 Ising spin model with the \/h x y/h vacancy superstructure for three representative 
sets of model parameters. The main results can be summarized as below. 

(i) The magnetic configuration of the groundstate in Case-I is the BSC state as observed 
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in realistic AyFe2-x^e2 compounds. While, the groundstate magnetic configuration in Case- 
II is the SAFM state as observed in other iron pnictides, though there are iron vacancies 
here. 

(ii) A dramatic enhancement in the transition temperature, Tc = 3.6238(4), is determined 
in Case-I. This value is almost three times as large as that in Case-II where = 1.2273. 



This result is quite compati 



the BSC {Tn ~ 550 K) 
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e with the corresponding magnetic transition temperatures of 



25j and SAFM {Tn ~ 100 - 200 K) 



10 



11| phases reported 



respectively in experiments. 

(iii) The nature of the magnetic transition in Case-I is revealed of the 2D Ising criticality. 
Significant deviation of the critical exponents from the 2D Ising criticality, reaching about 
10%, indicates that Case-II belongs to the Suzuki's weak universality class. 

(iv) Case-Ill is shown to be equivalent to the bathroom-tile Ising model which is exactly 
solvable. Good agreement between the numerical results and the exact solution demonstrates 
the validity of our numerical simulations on this class of complex systems. 

Finally, with a block-spin rotation and topological deformation, we show the absence and 
persistence of the magnetic frustration caused by the competing antiferromagnetic NNN 
coupling in Case-I and Case-II, respectively. The groundstate magnetic configuration for 
model parameters obtained from the first-principles calculations [28;i| falls to the same BSC 
state as in Case-I where the magnetic frustration is significantly suppressed as seen after 
the block-spin rotation. Combined with the numerical results, we conclude that the dra- 
matic enhancement of the transition temperature in the BSC state as observed in realistic 
materials y4yFe2-xSe2 should be mainly due to a combination effect of the perfect vacancy 
superstructure and the block lattice contraction. The latter in turn leads to the suppression 
of magnetic frustrations due to the ferromagnetic intrablock couplings and the dominating 
antiferromagnetic interblock NNN coupling. 
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TABLE I: The inverse transition temperatures and critical exponents obtained with the short-time 
dynamic approach are hsted for Case-I, Case-II, and Case-Ill, in comparison with those of 2D 
Ising model on the square lattice from literatures 
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52 



55|. Not all of the critical exponents are 



independent, and the scaling laws ^ /u + 213 /u = d and \ + 9z = d hold quite well within error bars 
in each case. 





Case-I 


Case-II 


Case-Ill 


2D Ising 


Ground-state Kc 


0.27595(3) 0.8148(1) 0.6952(1) 0.44069 


/3 


0.125(2) 


0.109(1) 


0.122(2) 


1/8 


V 


1.00(1) 


0.90(1) 


1.00(2) 


1 


z 


2.16(1) 


2.17(1) 


2.18(2) 


2.16(1) 




0.125(2) 


0.121(2) 


0.122(3) 


1/8 


Disordered 6 


0.186(2) 


0.167(1) 


0.189(1) 


0.191(1) 


7 


1.72(2) 


1.58(2) 


1.76(3) 


7/4 




1.74(2) 


1.75(2) 


1.77(2) 


7/4 


X 


1.59(1) 


1.65(1) 


1.60(1) 


1.59(1) 
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Block spin checkerboard Stripe-like 

FIG. 1: The magnetic structure of ^yFe2-a;Se2 from top-view is displayed. In the upper panel, 
the solid squares connecting the circles indicate the fundamental blocks with four Fe atoms at the 
corners. The proposed magnetic couplings (Ji^Ji) with solid lines and (J2, ^2) with dashed lines 
represent the NN and NNN couplings, respectively. In the lower panel, two magnetic configurations, 

the block spin checkerboard and stripe-like antifcrromagnetic states, are shown by open (Si = 1) 
and solid circles {Si = — 1). 
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(a) (b) 

FIG. 2: (a) The second moment of the pseudo-magnetization in the parallel tempering process is 
displayed for different sets of model parameters. Dashed line represents constant, M^'^\t) = 1. 
(b) Dynamic relaxation of the pseudo-magnetization is plotted for different temperatures. For 

clarity, the curves at Kc with different lattice sizes are shifted down. Dashed lines indicate power- 
law fits. 
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FIG. 3: (a) The logarithmic derivative of the pseudo-magnetization M{t,T) is displayed at Kc for 
Case-I (stars) and II (open squares). Dashed lines represent power-law fits, and solid lines indicate 
the fits with power-law correction. (b) The Binder cumulant U{t) is plotted with solid lines on a 
double- log scale for different lattice size L = 250, 500 and 1000 (from above). According to Eq. ([9]), 
data collapse is demonstrated at a fixed lattice size L = 500. Open circles and triangles correspond 
to -L = 250 and 1000, respectively. 
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FIG. 4: (a) The spatial correlation function C{r,t) is displayed on a log-log scale. Data collapse is 
demonstrated at a fixed t = 20480 MCS. Open circles, open triangles, stars, solid circles, and solid 
squares correspond to t = 20, 80, 320, 1280, and 5120, respectively. In the inset, the scaling function 
C{s)s^''^^ is shown on a linear-log scale. For clarity, the curve of Case-II is shifted right. (b) 
Dynamic relaxation of xi^) S{t) are plotted with solid lines for Case-I. Dashed lines represent 




(a) (b) 

FIG. 5: (a) The scaling function A{t/t') against is displayed for Case-I and II on a double- 

log scale. According to Eq. (llOp . data collapse is observed for different waiting time t' . (b) The 
time evolution of M{t) is plotted for Case-I and II with an initial magnetization tjiq = 0.01. The 
lattice size is L = 1000. In both (a) and (b), dashed lines show power-law fits. 
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(a) (b) 

FIG. 6: (a) M{t) in Case-Ill is plotted for different temperatures. For clarity, the curve at Kc 
is shifted down. Dashed line shows a power-law fit. (b) The square lattice with a \/5 x s/b 
vacancy superstructure and the bathroom-tile lattice are shown within the dashed squares. A 
same topological structure is revealed. 
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